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Abstract 

This paper reports investigations on the computation of mate- 
rial fronts in multi-fluid models using a Lagrange-Projection ap- 
proach. Various forms of the Projection step are considered. Par- 
ticular attention is paid to minimization of conservation errors. 

1 Introduction 

It is well-known that standard conservative discretizations of gas dy- 
namics equations in Eulerian coordinates generally develop non physi- 
cal pressure oscillations near contact discontinuities, and more generally 
near material fronts in multi-component flows. Several cures based on 
a local non conservative modification have been proposed. Let us quote 
for instance the hybrid algorithm derived by Karni in [TU] and the Two- 
Flux method proposed by Abgrall and Karni in |2 for multi-fluid flows. 
See also [1], |3] and the references therein. 

We investigate here a Lagrange-Projection type method to get rid of 
pressure oscillations. The basic motivation lies in the fact that oscil- 
lations do not exist in Lagrangian computations. It is then possible 
to clearly determine which operation in the projection step sparks off 
pressure oscillations. As in [TU], [2, a non conservative correction is pro- 
posed. It is based on a local pressure averaging and a random sampling 
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strategy on the mass fraction in order to strictly preserve isolated mate- 
rial fronts and get a statistical conservation property. Numerical results 
are proposed and compared with the Two-Flux method [2j. 



2 The model under consideration 

We consider a nonlinear partial differential equations model governing 
the flow of two species Si and E2 separated by a material interface. For 
instance, we focus on two perfect gases and we set Pi{pi, e^) — (7— l)piei, 
7i = Cp^i/Cy^i and T^{p^,e^) = ei/Cv,i where p^, pi, e^, 7i > 1, Ti, Cp^i > 
0, Cv^i > respectively denote the pressure, the density, the internal 
energy, the adiabatic coefficient, the temperature and the specific heats 
of Si, i = 1,2. The mixture density is given hy p = pi+ p2 and we adopt 
a Dalton's law for the mixture pressure p — pi(pi,ei) + P2{p2,e2)- We 
assume in addition that the two species evolve according to the same 
velocity u and are at thermal equilibrium, that is T = Ti(pi,ei) = 
12(^2,62). The mixture internal and total energies e and E are defined 
by pe = pici + ^2^2 and pE — \pv? + pe. Then, introducing the mass 

p = (7 - i)pe with 7 ^ 7(y) = ^ ' y — > i. 



fraction Y — pi/p, straightforward manipulations yield pe — pCyT with 

^ c„(r) = + (1 - y)a,2 and 



(2.1) 



yCp,i + (1 - Y)Cp 

YC^.l + (1 - Y)Cy,2 

In one-space dimension, the model under consideration writes 

' dtp + dx{pu) = 0, 
dt{pu) dx{pu^ +p) = 0, 
dt{pE) + d^{pE+p)u^Q, 
ydtpY + d^{pYu) = Q, 

and for the sake of conciseness, we set 

dtn{x,t) + dJ{u{x,t))=Q. (2.2) 

The flux function f flnds a natural definition with respect to the conser- 
vative unknowns u = (p, pu, pE, pY). Let us mention that (j2.2p is hyper- 
bohc with eigenvalues Ao(u) = u and A±(u) = u± c(u), c(u) — \fjplp, 
provided that p>0, 0<F<1 and p > 0. The characteristic field asso- 
ciated with Aq is linearly degenerate, leading to contact discontinuities 
or material fronts. The two extreme fields are genuinely nonlinear. 



3 Numerical schemes 

This section is devoted to the discretization of (|2.2I) . As already stated, 
a specific attention must be paid to the contact discontinuities to avoid 
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pressure oscihations. With this in mind, we first revisit the "Two-Flux 
Method" proposed by Abgrall and Karni [2 and then present a new 
numerical procedure based on a Lagrangian approach and a random 
sampling strategy. Comparisons will be proposed in section |4l 
Let us introduce a time step At > and a space step Aa; > that we 
assume to be constant for simplicity. We set A = At/ Ax and define the 
mesh interfaces ~ jAx for j G Z, and the intermediate times 

i" — nAt for n € N. In the sequel, u" denotes the approximate value of 
u at time t" and on the cell Cj — [xj_i/2, 2^^+1/2 [■ For n = and j € Z, 
we set u" = ^x'^l'/l uo(a;)(ix where Uo(a;) is the initial condition. 

3.1 The Two-Flux method revisited 

Aim of this section is to review the Two-Flux method proposed by Ab- 
grall and Karni [5]. Let us first recall that pressure oscillations do not 
systematically appear in single-fluid computations. Abgrall and Karni 
[5] then propose to replace any conservative multi-fluid strategy by a 
non conservative approach based on the definition of two single- fluid nu- 
merical fiuxes at each interface. We first recall the algorithm in details 
and then suggest a slight modification in order to lessen the conservation 
errors. This strategy will be used as a reference to assess the validity of 
the Lagrangian strategies proposed in the next subsection. 

The original algorithm. Let us consider a two-point numerical flux 
function g consistent with f . The Two-Flux method proposes to update 
the sequence (u")jgz into two steps, under an usual 1/2 CFL restriction. 

First step : evolution of p, pu and p {t" — )■ t"^^^) 

Let us first define v = {p, pu,p, Y) and the one-to-one mapping u — >■ v = 
v(u) thanks to the thermodynamics closures. Two interfacial numerical 
fluxes gj+i/2.L and gj+i/2,_R are then defined by 



where u^^ij, = with v^+i = {p]^^, {pu)'}^-^,p]^-^,YJ'), and 

u^_fl = u(v"_^) with v"^^ = {p'j,{pu)^,p'j,Yp^;^). In some sense, the 



mass fraction Y is then assumed to be the same on both side of each 
interface since Y — Y" is used for the computation of gj+i/2,L and 
Y = YJl^ for gj+i/2,fl. At last, we use gj+i/2,L, respectively gj+i/2,R, 
to update the conservative unknowns p, pu and pE on the cell j, resp. 
(j + 1). With clear notations, we get for all j G Z 



gj+i/2.L = g(u",uj 




(3.1) 




n+l- 




(3.2) 



ipE); 



)■ 
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Let us note that pi = pY is not concerned with p.2p . We simply set 
(pi)"^ = p^j^^^ X (Y is naturally kept constant in this step) and 
thenp;^+i- -p(uj+i-) with u]+'- - {p,pu,pE,p^)]+'-. 

Second step : evolution of pi {t""^^^ 

In this step, pi is evolved in a conservative way using the numerical flux 
function g while the values of p, pu and p are kept unchanged. Again 
with clear notations, the vector vj+i = y/+i) is 

then defined - p']+^- , {pu)]+^ = (P")^'", = pI^'^ 

y^"+i = {pi)]+'/p]+\ where 

(Pi)r-(/'i);-^(<i/2-gfi/2)- (3-3) 
with = g(u;',u;'+i). At last, we set u]+^ = u{v]+^). 

It is worth noticing that the Two-Flux method is not conservative on p 
and pu since the fluxes gj+i/2,L and gj+i/2,R defined at each interface 
are different as soon as Yj" 7^ ^j+i- K is not conservative on pE either, 
but by construction the Two-Flux method is conservative on pi. 

The associated quasi-conservative algorithm. It is actually clear 
from [To], [2] and the references therein that in standard conservative 
discretizations of p.2p . only the update formula of the total energy pE 
is responsible for the pressure oscillations. We are then tempted to pro- 
pose a quasi-conservative variant of the Two-Flux method such that only 
the total energy is treated in a non conservative form. For all j G Z, we 
simply replace p.2p by 

{pu);+'- = {pur; - - izA) 
(pi?);+^- = (pi?);-A(g;^/,,,-gf,/,,^), 

where again gj+1/2 ~ g(u", u"_^j). The second step is unchanged. 

3.2 The Lagrange-Projection approach with random 
samphng 

In this section, we propose a Lagrangian approach for approximating the 
solutions of (|2.2p . The general idea is to first solve this system in La- 
grangian coordinates, and then to come back to an Eulerian description 
of the flow with a projection step. Under its classical conservative form, 
the Lagrange-Projection method generates spurious oscillations near the 
material fronts. In order to remove these oscillations, we propose to 
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adapt the projection step (only). We begin with a description of the 
Lagrangian step and then recall, for the sake of clarity, the usual con- 
servative projection step (see for instance [S]). Again, an usual 1/2 CFL 
restriction is used. 



The Lagrangian step (t" 

In this step, ()2.2|) in written in Lagrangian coordinates and solved by an 
acoustic scheme (see for instance [7_), which gives 



where the velocity and the pressure at interfaces are defined by 



(3.5) 



1 (pc)" 



'j+1/2 
1/2, 



(3.6) 



The proposed local approximation ipc)"_^^-^^2 '^^ Lagrangian sound 
speed is (pc)"_^-^^2 = ^^^iiP'^)J Ap^)^+i) but other definitions may be 
found for instance in In this step, the grid points move at 

velocity 7/^_^i/2 so that p';+'- = l/r;+'', {puT+^- - p';+^- x wj+i-, 

{pE)^+^- = X and (pi)"+'" = p"+'" X define 

approximate values of u on a Lagrangian grid with mesh interfaces 

^i+1/2 = ^i+1/2 + '«"+l/2^*- 

The usual projection step — >■ t""^^) 

Aim of this step is to project the solution obtained at the end of the first 
step on the Eulerian grid defined by the mesh interfaces Xj+if2- Usu- 
ally, the choice is made to project the conservative vector u in order to 
obtain a conservative Lagrange-Projection scheme (see again [8]). More 
precisely, such a choice writes 

^"+1 = / ip''+^-{x)dx with ifi = p,pu,pE,pi = pY, (3.7) 



^ Ax 



\ -l/2if u", , ,„ > 0' 

or, with Ax* = x*^^/2 - and e{j, n) ^ i u^^^^ < 0' 



~ A^^ ^'^J' ^H"j + l/2'^J + l/2+£a,n) "j-l/2'Pj-l/2+eO-l,n)ij-- 
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What is wrong with this scheme ? The Lagrange-Projection ap- 
proach allows to precisely reveal the operation that makes the material 
fronts necessarily damaged by pressure oscillations. Let us indeed con- 
sider an isolated material front with uniform velocity and pressure pro- 
files : u{x,0) = u° > and p(a;,0) = while Y{x,0) = 1 if a; < 
and Y{x, 0) = if a; > 0. The density is also set to be uniform for 
simplicity : p{x,0) — . We first observe that this profile is clearly 
preserved in the Lagrangian step since by p.6p we have u^_|_i/2 ^ '^^ ^'^"^ 
?'i-l-i/2 ~ P° P-6P - Then, the projection procedure p.7p gives — p° 
and (/cm)] = so that m] — {pu)j/pj — After the first time 

iteration, the velocity profile is then still free of spurious oscillations. At 
last, using the property that this velocity is constant and positive, and 
focusing for instance on the cell of index j = 1, (|3.7p gives for ip = pE, pi 

(pE)\ = {pEt-u^\{{pEt-{pE)l) and = ^°A. 

The pressure is then given after easy calculations by 

72 ^ J- 71"-'- 

At this stage, there is no reason for p\ to equal p'^ From the very first 
time iteration, a pressure oscillation is then created. As an immediate 
consequence, the velocity profile will not remain uniform in the next time 
iteration and the numerical solution is damaged for good. 

It is then clear that the way the pressure is updated in the projection 
step (only) is responsible for the spurious oscillations in an usual conser- 
vative Lagrange-Projection scheme. We propose to modify this step. As 
for the Two-Flux method, the idea is to give up the conservation prop- 
erty in order to maintain uniform the pressure p across material fronts. 
Let us emphasize that the Lagrangian step is unchanged. 

The quasi-conservative p-projection step (t"^^^ — >■ t"^^ ) 
First of all, p, pu and pi still evolve according to p.7p so that the al- 
gorithm remains conservative on these variables. We will keep on using 
p.7p for — pE only for j not in a subset of Z defined below. On 
the contrary, the pressure p (instead of pE) is averaged for j g Z^: 

pT = -JT / p''^^-{x)dx. (3.8) 

Fori e Z;,wethensetv;'+i - {p']+\{puT+\p']+\Y;'+^) wilhY^+^ = 
{Pi)T'IpT\ and <+i = u(v"+i). 
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Definition of Zp. Up to our knowledge, the idea of averaging the 
pressure p in a Lagrange-Projection strategy first appeared in 5^. This 
way to proceed is clearly sufficient to remove the pressure oscillations 
near the material fronts if = Z. However, averaging the pressure 
p instead of the total energy pE for all j e Z gives a non conserva- 
tive scheme that is expected to provide discontinuous solutions violating 
the Rankine-Hugoniot conditions, see for instance Hou and LeFloch [9] 
(note however that here, the Lagrangian system associated with (|2.ip 
is actually treated in a conservative form, while the pressure averaging 
takes place in the projection step only). This was confirmed in prac- 
tice when considering solutions involving shocks with large amplitude. 
Then, in order to lessen the conservation errors, we propose to localize 
the non conservative treatment around the contact discontinuities set- 
ting Z; = [j e Z,max(|y/ - YJl^l \Y^^^ - > e} for a given e > 0. 
Following Karni [10], we will use e = 0.05 in practice. 

The quasi- conservative p-projection step with sampling (t"'~^^~ — >■ t"'~^^) 
The quasi-conservative p-projection step will be seen in the next section 
to properly compute large amplitude shock propagations. Localizing the 
averaging process of p nevertheless prevents the method from keeping 
strictly uniform the velocity and pressure profiles of an isolated material 
front, see Test A below. Indeed, note that since Z^ is generally a strict 
subset of Z due to the numerical diffusion on Y {i.e. Z^ C Z), an usual 
conservative treatment is still used on pE as soon as j not in Z^. This is 
sufficient to create pressure oscillations. In order to cure this problem, 
we propose to get rid of the numerical diffusion on Y so as to enforce the 
non conservative treatment (|3.8|) across an isolated material front. This 
objective is achieved when replacing the conservative updating formula 
(|3.7p for pi with random sampling strategy applied to Y (see also [4] 
and [S] for similar ideas) . More precisely, we consider an equidistributed 
random sequence (a,i) in (0, 1) (following Collela [6], we take in practice 
the celebrated van der Corput sequence), define x* — Xj_i/2 + On+iAt 
for all j 6 Z and set 







if 














if 


X* 




i 


3- 






if 


^1 



< X 



i-1/2' 
— — ^j+i/2' 



> X 



j+1/2- 



Then, we set {pi)"^^ — Pj^^^"^^ ^'^ ^^^^ the conservation of pi now 
holds only statistically. 



8 



C. Chalons, F. Coquel 



4 Numerical results 

We propose two numerical experiments with 71 ~ 1.4 and 72 ~ 1.6 
associated with a Riemann initial data. The left and right vectors v are 
denoted and and the initial discontinuity is at a: = 0.5. In the first 
simulation (Test A), wc consider the propagation of an isolated material 
interface with = (1, 1, 1, 1) and v^; — (0.1, 0.1, 1, 0). We take Ax — 
0.005 and plot the solutions at time t — 0.15. The second simulation 
(Test B) develops a strong shock due to a large initial pressure ratio. 
More precisely, we choose vl = (1,0,500, 1) and vj^ = (1,0,0.2,0). We 
take Ax = 0.00125 and plot the solutions at time t — 0.008. 
We observe that the Two-Flux method and the Lagrangian methods 
are in agreement with the exact solutions and give similar results. As 
expected, note that the Lagrangian approach without sampling does 
not strictly maintain uniform the pressure profile for Test A. Note also 
that the mass fraction Y is sharp when a random sampling is used. At 
last, the relative conservation error on pE (see for instance [5] for more 
details) for the Lagrangian approach with random sampling is actually 
less important and swings around 0.2% only. 




Figure 4.1: p (Test A) Figure 4.2: p (Test A) 



5 Concluding remarks 

We have investigated a Lagrange-Projection approach for computing ma- 
terial fronts in multi-fiuid models. We get similar results to the Two-Flux 
method 2_ with less important conservation errors on pE. Let us men- 
tion that other strategies, like for instance the one consisting in a local 
random sampling of u (instead of Y only) in the Lagrangian step, have 
been investigated. The results are not reported here. 

Acknowledgements. The authors are grateful for helpful discussions 



Computing material fronts with a Lagrange-Projection approach 9 




Figure 4.3: p (Test B) Figure 4.4: Y (Test B) 




Figure 4.5: Conservation errors (A)Figure 4.6: Conservation errors (B) 
and exchanges with P. Helluy and F. Lagoutiere. 
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